{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib ipympl"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "# The double pendulum problem\n",
        "\n",
        "This animation illustrates the double pendulum problem.\n",
        "\n",
        "Double pendulum formula translated from the C code at\n",
        "http://www.physics.usyd.edu.au/~wheat/dpend_html/solve_dpend.c\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "application/vnd.jupyter.widget-view+json": {
              "model_id": "912300e74e8543fca0a4b480b95a9395",
              "version_major": 2,
              "version_minor": 0
            },
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfQAAAGQCAYAAABYs5LGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAkVElEQVR4nO3df1TW9f3/8QfoxYWsIYnIBaWGWtIP8wdOB9tJmBhUp8nquCzP/HFMy+VOjlZpZ9Nha8zmstVs1qmkneVqnZPuR80OkdhZERbJyZyxNJNlXthwQErhNXl9/+jr9YkJBMrF+3q9ut/Ouc7p/eb9gsfz+uGj6+J9ccUYY4wAAIDVYr0OAAAAzhyFDgCAAyh0AAAcQKEDAOAACh0AAAdQ6AAAOIBCBwDAARQ6AAAOoNABAHAAhQ4AgAModAAAHEChAwDgAAodAAAHUOgAADiAQgcAwAEUOgAADqDQAQBwAIUOAIADKHQAABxAoQMA4AAKHQAAB1DoAAA4gEIHAMABFDoAAA6g0AEAcACFDgCAAyh0AAAcQKEDAOAACh0AAAdQ6AAAOIBCBwDAARQ6AAAOoNABAHAAhQ4AgAModAAAHEChAwDgAAodAAAHUOgAADiAQgcAwAEUOgAADqDQAQBwAIUOAIADKHQAABxAoQMA4AAKHQAAB1DoAAA4gEIHAMABFDoAAA6g0AEAcACFDgCAAyh0AAAc4Hmhv/zyy7r66quVnp6umJgYbdmy5QvXVFZWatKkSfL7/RozZozKysoinhMAgGjmeaEfO3ZM48eP1/r163t0/P79+3XVVVcpLy9PtbW1WrZsmW688Ua98MILEU4KAED0ijHGGK9DnBQTE6PNmzerqKioy2PuvPNOPffcc3r77bfD+2bPnq2mpiZt3bq1H1ICABB9PH+G3ltVVVXKz8/vsK+goEBVVVUeJQIAwHsDvQ7QW8FgUKmpqR32paamqqWlRZ988okGDRrU6bq2tja1tbWFt9vb23XkyBElJycrJiYmopkBANHJGKOPP/5Y6enpio217jluB9YV+ukqLS1VSUmJ1zEAAFHoX//6l84991yvY5wR6wo9EAiooaGhw76GhgYlJiZ2+excklasWKHi4uLwdnNzs0aMGKF//vOfGjJkSMTyRkooFNK2bduUl5cnn8/ndZxesz2/ZP8MtueX7J/B9vyS/TMcOXJEF1xwgb761a96HeWMWVfo2dnZev755zvsKy8vV3Z2drfr/H6//H7/KfuHDBmi5OTkPs3YH0KhkBISEpScnGzlg8j2/JL9M9ieX7J/BtvzS27MIMmJX716/guDo0ePqra2VrW1tZI+e1tabW2t6uvrJX32zHru3Lnh42+++Wa99957uuOOO/TOO+/ooYce0h//+Ef98Ic/9CI+AABRwfNCf+ONNzRx4kRNnDhRklRcXKyJEydq5cqVkqRDhw6Fy12SMjIy9Nxzz6m8vFzjx4/Xr371Kz366KMqKCjwJD8AANHA85fcc3Nz1d1b4Tv7K3C5ubnauXNnBFMBAGAXz5+hAwCAM0ehAwDgAAodAAAHUOgAADiAQgcAwAEUOgAADqDQAQBwAIUOAIADKHQAABxAoQMA4AAKHQAAB1DoAAA4gEIHAMABFDoAAA6g0AEAcACFDgCAAyh0AAAcQKEDAOAACh0AAAdQ6AAAOIBCBwDAARQ6AAAOoNABAHAAhQ4AgAModAAAHEChAwDgAAodAAAHUOgAADiAQgcAwAEUOgAADqDQAQBwAIUOAIADKHQAABxAoQMA4AAKHQAAB1DoAAA4gEIHAMABFDoAAA6ImkJfv369zjvvPMXHx2vq1KnasWNHl8eWlZUpJiamwyU+Pr4f0wIAEF2iotCffvppFRcXa9WqVXrzzTc1fvx4FRQU6PDhw12uSUxM1KFDh8KXAwcO9GNiAACiS1QU+n333adFixZpwYIFuuiii7RhwwYlJCTo8ccf73JNTEyMAoFA+JKamtqPiQEAiC4DvQ5w/Phx1dTUaMWKFeF9sbGxys/PV1VVVZfrjh49qpEjR6q9vV2TJk3Sz3/+c1188cVdHt/W1qa2trbwdktLiyQpFAopFAr1wST962RmG7NL9ueX7J/B9vyS/TPYnl+yfwZbc3cmxhhjvAzw4Ycf6pxzztGrr76q7Ozs8P477rhD27dvV3V19Slrqqqq9O677+rSSy9Vc3Oz1q5dq5dfflm7d+/Wueee2+nP+elPf6qSkpJT9m/atEkJCQl9NxAAwBqtra264YYb1NzcrMTERK/jnBHPn6Gfjuzs7A7ln5OTowsvvFAPP/yw7r777k7XrFixQsXFxeHtlpYWDR8+XHl5eUpOTo545r4WCoVUXl6uGTNmyOfzeR2n12zPL9k/g+35JftnsD2/ZP8MjY2NXkfoM54X+tChQzVgwAA1NDR02N/Q0KBAINCj7+Hz+TRx4kTt3bu3y2P8fr/8fn+na228E55Efu/ZPoPt+SX7Z7A9v2TvDDZm7ornJ8XFxcUpKytLFRUV4X3t7e2qqKjo8Cy8OydOnNCuXbuUlpYWqZgAAEQ1z5+hS1JxcbHmzZunyZMna8qUKbr//vt17NgxLViwQJI0d+5cnXPOOSotLZUkrV69Wl//+tc1ZswYNTU16Ze//KUOHDigG2+80csxAADwTFQU+nXXXaePPvpIK1euVDAY1IQJE7R169bwW9Hq6+sVG/t/Lyb85z//0aJFixQMBnX22WcrKytLr776qi666CKvRgAAwFNRUeiStHTpUi1durTTr1VWVnbYXrdundatW9cPqQAAsIPnv0MHAABnjkIHAMABFDoAAA6g0AEAcACFDgCAAyh0AAAcQKEDAOAACh0AAAdQ6AAAOIBCBwDAARQ6AAAOoNABAHAAhQ4AgAModAAAHEChAwDgAAodAAAHUOgAADiAQgcAwAEUOgAADqDQAQBwAIUOAIADKHQAABxAoQMA4AAKHQAAB1DoAAA4gEIHAMABFDoAAA6g0AEAcACFDgCAAyh0AAAcQKEDAOAACh0AAAdQ6AAAOIBCBwDAARQ6AAAOoNABAHAAhQ4AgAOiptDXr1+v8847T/Hx8Zo6dap27NjR7fHPPPOMMjMzFR8fr3Hjxun555/vp6QAAESfqCj0p59+WsXFxVq1apXefPNNjR8/XgUFBTp8+HCnx7/66qu6/vrrtXDhQu3cuVNFRUUqKirS22+/3c/JAQCIDlFR6Pfdd58WLVqkBQsW6KKLLtKGDRuUkJCgxx9/vNPjf/3rX6uwsFC33367LrzwQt19992aNGmSfvOb3/RzcgAAooPnhX78+HHV1NQoPz8/vC82Nlb5+fmqqqrqdE1VVVWH4yWpoKCgy+MBAHDdQK8D/Pvf/9aJEyeUmpraYX9qaqreeeedTtcEg8FOjw8Gg13+nLa2NrW1tYW3W1paJEmhUEihUOh043vmZGYbs0v255fsn8H2/JL9M9ieX7J/Bltzd8bzQu8vpaWlKikpOWX/tm3blJCQ4EGivlFeXu51hDNie37J/hlszy/ZP4Pt+SV7Z2htbfU6Qp/xvNCHDh2qAQMGqKGhocP+hoYGBQKBTtcEAoFeHS9JK1asUHFxcXi7paVFw4cPV15enpKTk89gAm+EQiGVl5drxowZ8vl8XsfpNdvzS/bPYHt+yf4ZbM8v2T9DY2Oj1xH6jOeFHhcXp6ysLFVUVKioqEiS1N7eroqKCi1durTTNdnZ2aqoqNCyZcvC+8rLy5Wdnd3lz/H7/fL7/afs9/l8Vt4JTyK/92yfwfb8kv0z2J5fsncGGzN3xfNCl6Ti4mLNmzdPkydP1pQpU3T//ffr2LFjWrBggSRp7ty5Ouecc1RaWipJuvXWWzVt2jT96le/0lVXXaWnnnpKb7zxhh555BEvxwAAwDNRUejXXXedPvroI61cuVLBYFATJkzQ1q1bwye+1dfXKzb2/07Iz8nJ0aZNm/TjH/9Yd911l84//3xt2bJFl1xyiVcjAADgqagodElaunRply+xV1ZWnrJv1qxZmjVrVoRTAQBgB8/fhw4AAM4chQ4AgAModAAAHEChAwDgAAodAAAHUOgAADiAQgcAwAEUOgAADqDQAQBwAIUOAIADKHQAABxAoQMA4AAKHQAAB1DoAAA4gEIHAMABFDoAAA6g0AEAcACFDgCAAyh0AAAcQKEDAOAACh0AAAdQ6AAAOIBCBwDAARQ6AAAOoNABAHAAhQ4AgAModAAAHEChAwDgAAodAAAHUOgAADiAQgcAwAEUOgAADqDQAQBwAIUOAIADKHQAABxAoQMA4AAKHQAAB1DoAAA4wPNCP3LkiObMmaPExEQlJSVp4cKFOnr0aLdrcnNzFRMT0+Fy880391NiAACiz0CvA8yZM0eHDh1SeXm5QqGQFixYoMWLF2vTpk3drlu0aJFWr14d3k5ISIh0VAAAopanhb5nzx5t3bpVr7/+uiZPnixJevDBB3XllVdq7dq1Sk9P73JtQkKCAoFAf0UFACCqeVroVVVVSkpKCpe5JOXn5ys2NlbV1dX6zne+0+XaJ598Ur///e8VCAR09dVX6yc/+Um3z9Lb2trU1tYW3m5paZEkhUIhhUKhPpimf53MbGN2yf78kv0z2J5fsn8G2/NL9s9ga+7OeFrowWBQw4YN67Bv4MCBGjJkiILBYJfrbrjhBo0cOVLp6el66623dOedd6qurk7PPvtsl2tKS0tVUlJyyv5t27ZZ/XJ9eXm51xHOiO35JftnsD2/ZP8MtueX7J2htbXV6wh9JiKFvnz5cq1Zs6bbY/bs2XPa33/x4sXh/x43bpzS0tI0ffp07du3T6NHj+50zYoVK1RcXBzebmlp0fDhw5WXl6fk5OTTzuKVUCik8vJyzZgxQz6fz+s4vWZ7fsn+GWzPL9k/g+35JftnaGxs9DpCn4lIod92222aP39+t8eMGjVKgUBAhw8f7rD/v//9r44cOdKr349PnTpVkrR3794uC93v98vv95+y3+fzWXknPIn83rN9BtvzS/bPYHt+yd4ZbMzclYgUekpKilJSUr7wuOzsbDU1NammpkZZWVmSpJdeeknt7e3hku6J2tpaSVJaWtpp5QUAwHaevg/9wgsvVGFhoRYtWqQdO3bolVde0dKlSzV79uzwGe4HDx5UZmamduzYIUnat2+f7r77btXU1Oj999/Xn//8Z82dO1eXXXaZLr30Ui/HAQDAM57/YZknn3xSmZmZmj59uq688kp985vf1COPPBL+eigUUl1dXfjEhbi4OL344ou6/PLLlZmZqdtuu03XXnut/vKXv3g1AgAAnvP8D8sMGTKk2z8ic95558kYE94ePny4tm/f3h/RAACwhufP0AEAwJmj0AEAcACFDgCAAyh0AAAcQKEDAOAACh0AAAdQ6AAAOIBCBwDAARQ6AAAOoNABAHAAhQ4AgAModAAAHEChAwDgAAodAAAHUOgAADiAQgcAwAEUOgAADqDQAQBwAIUOAIADKHQAABxAoQMA4AAKHQAAB1DoAAA4gEIHAMABFDoAAA6g0AEAcACFDgCAAyh0AAAcQKEDAOAACh0AAAdQ6AAAOIBCBwDAARQ6AAAOoNABAHAAhQ4AgAModAAAHEChAwDgAM8L/Z577lFOTo4SEhKUlJTUozXGGK1cuVJpaWkaNGiQ8vPz9e6770Y2KAAAUczzQj9+/LhmzZqlJUuW9HjNvffeqwceeEAbNmxQdXW1vvKVr6igoECffvppBJMCABC9BnodoKSkRJJUVlbWo+ONMbr//vv14x//WDNnzpQk/e53v1Nqaqq2bNmi2bNnRyoqAABRy/NC7639+/crGAwqPz8/vG/w4MGaOnWqqqqquiz0trY2tbW1hbdbWlokSaFQSKFQKLKhI+BkZhuzS/bnl+yfwfb8kv0z2J5fsn8GW3N3xrpCDwaDkqTU1NQO+1NTU8Nf60xpaWn41YDP27ZtmxISEvo2ZD8qLy/3OsIZsT2/ZP8MtueX7J/B9vySvTO0trZ6HaHPRKTQly9frjVr1nR7zJ49e5SZmRmJH9+pFStWqLi4OLzd0tKi4cOHKy8vT8nJyf2Wo6+EQiGVl5drxowZ8vl8XsfpNdvzS/bPYHt+yf4ZbM8v2T9DY2Oj1xH6TEQK/bbbbtP8+fO7PWbUqFGn9b0DgYAkqaGhQWlpaeH9DQ0NmjBhQpfr/H6//H7/Kft9Pp+Vd8KTyO8922ewPb9k/wy255fsncHGzF2JSKGnpKQoJSUlEt9aGRkZCgQCqqioCBd4S0uLqqure3WmPAAALvH8bWv19fWqra1VfX29Tpw4odraWtXW1uro0aPhYzIzM7V582ZJUkxMjJYtW6af/exn+vOf/6xdu3Zp7ty5Sk9PV1FRkUdTAADgLc9Pilu5cqWeeOKJ8PbEiRMlfXayWm5uriSprq5Ozc3N4WPuuOMOHTt2TIsXL1ZTU5O++c1vauvWrYqPj+/X7AAARAvPC72srOwL34NujOmwHRMTo9WrV2v16tURTAYAgD08f8kdAACcOQodAAAHUOgAADiAQgcAwAEUOgAADqDQAQBwAIUOAIADKHQAABxAoQMA4AAKHQAAB1DoAAA4gEIHAMABFDoAAA6g0AEAcACFDgCAAyh0AAAcQKEDAOAACh0AAAdQ6AAAOIBCBwDAARQ6AAAOoNABAHAAhQ4AgAModAAAHEChAwDgAAodAAAHUOgAADiAQgcAwAEUOgAADqDQAQBwAIUOAIADKHQAABxAoQMA4AAKHQAAB1DoAAA4gEIHAMABFDoAAA7wvNDvuece5eTkKCEhQUlJST1aM3/+fMXExHS4FBYWRjYoAABRbKDXAY4fP65Zs2YpOztbjz32WI/XFRYWauPGjeFtv98fiXgAAFjB80IvKSmRJJWVlfVqnd/vVyAQiEAiAADs4/lL7qersrJSw4YN09ixY7VkyRI1NjZ6HQkAAM94/gz9dBQWFuqaa65RRkaG9u3bp7vuuktXXHGFqqqqNGDAgE7XtLW1qa2tLbzd0tIiSQqFQgqFQv2Suy+dzGxjdsn+/JL9M9ieX7J/BtvzS/bPYGvuzsQYY0xff9Ply5drzZo13R6zZ88eZWZmhrfLysq0bNkyNTU19frnvffeexo9erRefPFFTZ8+vdNjfvrTn4Zf3v+8TZs2KSEhodc/EwBgv9bWVt1www1qbm5WYmKi13HOSEQK/aOPPvrCl8BHjRqluLi48PaZFLokpaSk6Gc/+5luuummTr/e2TP04cOH69ChQ0pOTj6tn+mlUCik8vJyzZgxQz6fz+s4vWZ7fsn+GWzPL9k/g+35JftnaGxsVFpamhOFHpGX3FNSUpSSkhKJb92pDz74IHyjdMXv93d6JrzP57PyTngS+b1n+wy255fsn8H2/JK9M9iYuSuenxRXX1+v2tpa1dfX68SJE6qtrVVtba2OHj0aPiYzM1ObN2+WJB09elS33367XnvtNb3//vuqqKjQzJkzNWbMGBUUFHg1BgAAnvL8pLiVK1fqiSeeCG9PnDhRkrRt2zbl5uZKkurq6tTc3CxJGjBggN566y098cQTampqUnp6ui6//HLdfffdvBcdAPCl5Xmhl5WVfeF70D//a/5BgwbphRdeiHAqAADs4vlL7gAA4MxR6AAAOIBCBwDAARQ6AAAOoNABAHAAhQ4AgAModAAAHEChAwDgAAodAAAHUOgAADiAQgcAwAEUOgAADqDQAQBwAIUOAIADKHQAABxAoQMA4AAKHQAAB1DoAAA4gEIHAMABFDoAAA6g0AEAcACFDgCAAyh0AAAcQKEDAOAACh0AAAdQ6AAAOIBCBwDAARQ6AAAOoNABAHAAhQ4AgAModAAAHEChAwDgAAodAAAHUOgAADiAQgcAwAEUOgAADqDQAQBwAIUOAIADPC30999/XwsXLlRGRoYGDRqk0aNHa9WqVTp+/Hi36z799FPdcsstSk5O1llnnaVrr71WDQ0N/ZQaAIDo42mhv/POO2pvb9fDDz+s3bt3a926ddqwYYPuuuuubtf98Ic/1F/+8hc988wz2r59uz788ENdc801/ZQaAIDoM9DLH15YWKjCwsLw9qhRo1RXV6ff/va3Wrt2badrmpub9dhjj2nTpk361re+JUnauHGjLrzwQr322mv6+te/3i/ZAQCIJp4Wemeam5s1ZMiQLr9eU1OjUCik/Pz88L7MzEyNGDFCVVVVXRZ6W1ub2traOvwcSTpy5EgfJe9foVBIra2tamxslM/n8zpOr9meX7J/BtvzS/bPYHt+yf4ZTnaAMcbjJGcuqgp97969evDBB7t8di5JwWBQcXFxSkpK6rA/NTVVwWCwy3WlpaUqKSk5Zf8FF1xw2nkBAG5obGzU4MGDvY5xRiJS6MuXL9eaNWu6PWbPnj3KzMwMbx88eFCFhYWaNWuWFi1a1OeZVqxYoeLi4vB2U1OTRo4cqfr6eitvxJaWFg0fPlz/+te/lJiY6HWcXrM9v2T/DLbnl+yfwfb8kv0zNDc3a8SIEd2+MmyLiBT6bbfdpvnz53d7zKhRo8L//eGHHyovL085OTl65JFHul0XCAR0/PhxNTU1dXiW3tDQoEAg0OU6v98vv99/yv7BgwdbeSc8KTExkfwes30G2/NL9s9ge37J/hliY+1/F3dECj0lJUUpKSk9OvbgwYPKy8tTVlaWNm7c+IVXalZWlnw+nyoqKnTttddKkurq6lRfX6/s7Owzzg4AgI08/V+SgwcPKjc3VyNGjNDatWv10UcfKRgMdvhd+MGDB5WZmakdO3ZI+uwZ9cKFC1VcXKxt27appqZGCxYsUHZ2Nme4AwC+tDw9Ka68vFx79+7V3r17de6553b42skzDkOhkOrq6tTa2hr+2rp16xQbG6trr71WbW1tKigo0EMPPdSrn+33+7Vq1apOX4a3Afm9Z/sMtueX7J/B9vyS/TPYnv/zYowL5+oDAPAlZ/9ZAAAAgEIHAMAFFDoAAA6g0AEAcMCXotBd+ZjWe+65Rzk5OUpISDjlT992Zf78+YqJielw+fwH4vSn08lvjNHKlSuVlpamQYMGKT8/X++++25kg3bjyJEjmjNnjhITE5WUlKSFCxfq6NGj3a7Jzc095Ta4+eab+yXv+vXrdd555yk+Pl5Tp04Nv/2zK88884wyMzMVHx+vcePG6fnnn++XnN3pzQxlZWWnXNfx8fH9mLajl19+WVdffbXS09MVExOjLVu2fOGayspKTZo0SX6/X2PGjFFZWVnEc3alt/krKytPuf5jYmK6/bPckVRaWqqvfe1r+upXv6phw4apqKhIdXV1X7guGh8HPfGlKHRXPqb1+PHjmjVrlpYsWdKrdYWFhTp06FD48oc//CFCCbt3OvnvvfdePfDAA9qwYYOqq6v1la98RQUFBfr0008jmLRrc+bM0e7du1VeXq6//vWvevnll7V48eIvXLdo0aIOt8G9994b8axPP/20iouLtWrVKr355psaP368CgoKdPjw4U6Pf/XVV3X99ddr4cKF2rlzp4qKilRUVKS333474lm70tsZpM/+Ytnnr+sDBw70Y+KOjh07pvHjx2v9+vU9On7//v266qqrlJeXp9raWi1btkw33nijXnjhhQgn7Vxv859UV1fX4TYYNmxYhBJ2b/v27brlllv02muvqby8XKFQSJdffrmOHTvW5ZpofBz0mPmSuvfee01GRkaXX29qajI+n88888wz4X179uwxkkxVVVV/ROzSxo0bzeDBg3t07Lx588zMmTMjmqe3epq/vb3dBAIB88tf/jK8r6mpyfj9fvOHP/whggk7949//MNIMq+//np439/+9jcTExNjDh482OW6adOmmVtvvbUfEnY0ZcoUc8stt4S3T5w4YdLT001paWmnx3/3u981V111VYd9U6dONTfddFNEc3antzP05rHR3ySZzZs3d3vMHXfcYS6++OIO+6677jpTUFAQwWQ905P827ZtM5LMf/7zn37J1FuHDx82ksz27du7PCYaHwc99aV4ht6ZM/2YVptUVlZq2LBhGjt2rJYsWaLGxkavI/XI/v37FQwGO9wGgwcP1tSpUz25DaqqqpSUlKTJkyeH9+Xn5ys2NlbV1dXdrn3yySc1dOhQXXLJJVqxYkWHP5QUCcePH1dNTU2H6y42Nlb5+fldXndVVVUdjpekgoICz+7vpzODJB09elQjR47U8OHDNXPmTO3evbs/4vaJaLsNTteECROUlpamGTNm6JVXXvE6TtjJj83u7t9+m2+DqPr41P4SyY9pjTaFhYW65pprlJGRoX379umuu+7SFVdcoaqqKg0YMMDreN06eT2npqZ22O/VbRAMBk956XDgwIEaMmRIt3luuOEGjRw5Uunp6Xrrrbd05513qq6uTs8++2zEsv773//WiRMnOr3u3nnnnU7XBIPBqLmupdObYezYsXr88cd16aWXqrm5WWvXrlVOTo527959yl+jjEZd3QYtLS365JNPNGjQII+S9UxaWpo2bNigyZMnq62tTY8++qhyc3NVXV2tSZMmeZqtvb1dy5Yt0ze+8Q1dcsklXR4XbY+D3rD6Gfry5cs7PQHj85f/feBH+mNae+t0ZuiN2bNn69vf/rbGjRunoqIi/fWvf9Xrr7+uyspKK/L3h0jPsHjxYhUUFGjcuHGaM2eOfve732nz5s3at29fH04BScrOztbcuXM1YcIETZs2Tc8++6xSUlL08MMPex3tS2Hs2LG66aablJWVpZycHD3++OPKycnRunXrvI6mW265RW+//baeeuopr6NEjNXP0KPxY1p7q7cznKlRo0Zp6NCh2rt3r6ZPn37G3y+S+U9ezw0NDUpLSwvvb2ho0IQJE07re3ampzMEAoFTTsb673//qyNHjvTqPjF16lRJn71SNHr06F7n7YmhQ4dqwIABp7wro7v7byAQ6NXxkXY6M/wvn8+niRMnau/evZGI2Oe6ug0SExOj/tl5V6ZMmaK///3vnmZYunRp+CTWL3qlJtoeB71hdaG78DGtvZmhL3zwwQdqbGzsUJBnIpL5MzIyFAgEVFFRES7wlpYWVVdX9/pM/+70dIbs7Gw1NTWppqZGWVlZkqSXXnpJ7e3t4ZLuidraWknqs9ugM3FxccrKylJFRYWKiookffaSY0VFhZYuXdrpmuzsbFVUVGjZsmXhfeXl5Z59LPHpzPC/Tpw4oV27dunKK6+MYNK+k52dfcpbpLy8DfpCbW1tRO/r3THG6Ac/+IE2b96syspKZWRkfOGaaHsc9IrXZ+X1hw8++MCMGTPGTJ8+3XzwwQfm0KFD4cvnjxk7dqyprq4O77v55pvNiBEjzEsvvWTeeOMNk52dbbKzs70YwRhjzIEDB8zOnTtNSUmJOeuss8zOnTvNzp07zccffxw+ZuzYsebZZ581xhjz8ccfmx/96EemqqrK7N+/37z44otm0qRJ5vzzzzeffvpp1Oc3xphf/OIXJikpyfzpT38yb731lpk5c6bJyMgwn3zySb/nN8aYwsJCM3HiRFNdXW3+/ve/m/PPP99cf/314a//7/1o7969ZvXq1eaNN94w+/fvN3/605/MqFGjzGWXXRbxrE899ZTx+/2mrKzM/OMf/zCLFy82SUlJJhgMGmOM+d73vmeWL18ePv6VV14xAwcONGvXrjV79uwxq1atMj6fz+zatSviWbvS2xlKSkrMCy+8YPbt22dqamrM7NmzTXx8vNm9e7cn+T/++OPw/VySue+++8zOnTvNgQMHjDHGLF++3Hzve98LH//ee++ZhIQEc/vtt5s9e/aY9evXmwEDBpitW7dakX/dunVmy5Yt5t133zW7du0yt956q4mNjTUvvviiJ/mXLFliBg8ebCorKzv8u9/a2ho+xobHQU99KQp948aNRlKnl5P2799vJJlt27aF933yySfm+9//vjn77LNNQkKC+c53vtPhfwL627x58zqd4fOZJZmNGzcaY4xpbW01l19+uUlJSTE+n8+MHDnSLFq0KPyPYbTnN+azt6795Cc/Mampqcbv95vp06eburq6/g///zU2Nprrr7/enHXWWSYxMdEsWLCgw/+Q/O/9qL6+3lx22WVmyJAhxu/3mzFjxpjbb7/dNDc390veBx980IwYMcLExcWZKVOmmNdeey38tWnTppl58+Z1OP6Pf/yjueCCC0xcXJy5+OKLzXPPPdcvObvTmxmWLVsWPjY1NdVceeWV5s033/Qg9WdOvo3rfy8nM8+bN89MmzbtlDUTJkwwcXFxZtSoUR0eD/2tt/nXrFljRo8ebeLj482QIUNMbm6ueemll7wJb0yX/+5//jq15XHQE3x8KgAADrD6LHcAAPAZCh0AAAdQ6AAAOIBCBwDAARQ6AAAOoNABAHAAhQ4AgAModAAAHEChAwDgAAodAAAHUOgAADiAQgcAwAEUOgAADqDQAQBwAIUOAIADKHQAABxAoQMA4AAKHQAAB1DoAAA4gEIHAMABFDoAAA6g0AEAcACFDgCAAyh0AAAcQKEDAOAACh0AAAdQ6AAAOIBCBwDAARQ6AAAOoNABAHAAhQ4AgAP+Hw8uFz78i1lMAAAAAElFTkSuQmCC",
            "text/html": [
              "\n",
              "            <div style=\"display: inline-block;\">\n",
              "                <div class=\"jupyter-widgets widget-label\" style=\"text-align: center;\">\n",
              "                    Figure\n",
              "                </div>\n",
              "                <img src='' width=500.0/>\n",
              "            </div>\n",
              "        "
            ],
            "text/plain": [
              "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "from numpy import sin, cos\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "import scipy.integrate as integrate\n",
        "import matplotlib.animation as animation\n",
        "from collections import deque\n",
        "\n",
        "G = 9.8  # acceleration due to gravity, in m/s^2\n",
        "L1 = 1.0  # length of pendulum 1 in m\n",
        "L2 = 1.0  # length of pendulum 2 in m\n",
        "L = L1 + L2  # maximal length of the combined pendulum\n",
        "M1 = 1.0  # mass of pendulum 1 in kg\n",
        "M2 = 1.0  # mass of pendulum 2 in kg\n",
        "t_stop = 5  # how many seconds to simulate\n",
        "history_len = 500  # how many trajectory points to display\n",
        "\n",
        "\n",
        "def derivs(state, t):\n",
        "\n",
        "    dydx = np.zeros_like(state)\n",
        "    dydx[0] = state[1]\n",
        "\n",
        "    delta = state[2] - state[0]\n",
        "    den1 = (M1 + M2) * L1 - M2 * L1 * cos(delta) * cos(delta)\n",
        "    dydx[1] = (\n",
        "        M2 * L1 * state[1] * state[1] * sin(delta) * cos(delta)\n",
        "        + M2 * G * sin(state[2]) * cos(delta)\n",
        "        + M2 * L2 * state[3] * state[3] * sin(delta)\n",
        "        - (M1 + M2) * G * sin(state[0])\n",
        "    ) / den1\n",
        "\n",
        "    dydx[2] = state[3]\n",
        "\n",
        "    den2 = (L2 / L1) * den1\n",
        "    dydx[3] = (\n",
        "        -M2 * L2 * state[3] * state[3] * sin(delta) * cos(delta)\n",
        "        + (M1 + M2) * G * sin(state[0]) * cos(delta)\n",
        "        - (M1 + M2) * L1 * state[1] * state[1] * sin(delta)\n",
        "        - (M1 + M2) * G * sin(state[2])\n",
        "    ) / den2\n",
        "\n",
        "    return dydx\n",
        "\n",
        "\n",
        "# create a time array from 0..t_stop sampled at 0.02 second steps\n",
        "dt = 0.02\n",
        "t = np.arange(0, t_stop, dt)\n",
        "\n",
        "# th1 and th2 are the initial angles (degrees)\n",
        "# w10 and w20 are the initial angular velocities (degrees per second)\n",
        "th1 = 120.0\n",
        "w1 = 0.0\n",
        "th2 = -10.0\n",
        "w2 = 0.0\n",
        "\n",
        "# initial state\n",
        "state = np.radians([th1, w1, th2, w2])\n",
        "\n",
        "# integrate your ODE using scipy.integrate.\n",
        "y = integrate.odeint(derivs, state, t)\n",
        "\n",
        "x1 = L1 * sin(y[:, 0])\n",
        "y1 = -L1 * cos(y[:, 0])\n",
        "\n",
        "x2 = L2 * sin(y[:, 2]) + x1\n",
        "y2 = -L2 * cos(y[:, 2]) + y1\n",
        "\n",
        "fig = plt.figure(figsize=(5, 4))\n",
        "ax = fig.add_subplot(autoscale_on=False, xlim=(-L, L), ylim=(-L, 1.0))\n",
        "ax.set_aspect(\"equal\")\n",
        "ax.grid()\n",
        "\n",
        "(line,) = ax.plot([], [], \"o-\", lw=2)\n",
        "(trace,) = ax.plot([], [], \".-\", lw=1, ms=2)\n",
        "time_template = \"time = %.1fs\"\n",
        "time_text = ax.text(0.05, 0.9, \"\", transform=ax.transAxes)\n",
        "history_x, history_y = deque(maxlen=history_len), deque(maxlen=history_len)\n",
        "\n",
        "\n",
        "def animate(i):\n",
        "    thisx = [0, x1[i], x2[i]]\n",
        "    thisy = [0, y1[i], y2[i]]\n",
        "\n",
        "    if i == 0:\n",
        "        history_x.clear()\n",
        "        history_y.clear()\n",
        "\n",
        "    history_x.appendleft(thisx[2])\n",
        "    history_y.appendleft(thisy[2])\n",
        "\n",
        "    line.set_data(thisx, thisy)\n",
        "    trace.set_data(history_x, history_y)\n",
        "    time_text.set_text(time_template % (i * dt))\n",
        "    return line, trace, time_text\n",
        "\n",
        "\n",
        "ani = animation.FuncAnimation(fig, animate, len(y), interval=dt * 1000, blit=True)\n",
        "plt.show()\n"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3.8.10 ('venv': venv)",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.8.10"
    },
    "vscode": {
      "interpreter": {
        "hash": "e82bdbec295f1d416caad911bfe8299822751c36b2e078ad810889afff869327"
      }
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}
